#include <gmp.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

/* It should be enough :D */
#define DEFAULT_N_ITERATIONS 25
#define DEFAULT_PRECISION    35000000

int main(int argc, char *argv[])
{
    int n_iter = DEFAULT_N_ITERATIONS;
    size_t prec = DEFAULT_PRECISION;

    /* First argument overrides the default number of iterations */
    if(argc > 1)
        sscanf(argv[1], "%d", &n_iter);

    /* Second argument overrides the default precision in bits */
    if(argc > 2)
        sscanf(argv[2], "%zu", &prec);

    /* Multi precision floating points variables */
    mpf_t pi;
    mpf_t a, b, t, p;
    mpf_t a0, b0, t0, p0;

    mpf_set_default_prec(prec);

    /* Initialize all variables */
    mpf_init(pi);

    mpf_init_set_d(a, 1.0);

    mpf_init_set_ui(b, 2);
    mpf_sqrt(b, b);
    mpf_div(b, a, b);

    mpf_init_set_d(t, 0.25);
    mpf_init_set_d(p, 1.0);

    mpf_init(a0);
    mpf_init(b0);
    mpf_init(t0);
    mpf_init(p0);

    int i;
    for(i = 1; i < n_iter; i++)
    {
        mpf_set(a0, a);
        mpf_set(b0, b);
        mpf_set(t0, t);
        mpf_set(p0, p);

        /* Calculate a(n+1) from a(n) and b(n) */
        /* a(n+1) = { a(n) + b(n) } / 2 */
        mpf_add(a, a0, b0);
        mpf_div_ui(a, a, 2);

        fprintf(stderr, "a(%d)\n", i);

        /* Calculate b(n+1) from a(n) and b(n) */
        /* b(n+1) = sqrt{ a(n) * b(n) } */
        mpf_mul(b, a0, b0);
        mpf_sqrt(b, b);

        fprintf(stderr, "b(%d)\n", i);

        /* Calculate t(n+1) from a(n+1), a(n), p(n) and t(n) */
        /* t(n+1) = p(n) * { a(n) - a(n+1) } ^ 2 - t(n) */
        mpf_sub(t, a0, a);
        mpf_mul(t, t, t);
        mpf_mul(t, t, p0);
        mpf_sub(t, t0, t);

        fprintf(stderr, "t(%d)\n", i);

        /* Calculate p(n+1) from p(n)*/
        /* p(n+1) = 2 * p(n) */
        mpf_mul_ui(p, p0, 2);

        fprintf(stderr, "p(%d)\n", i);
    }

    /* Calculate pi */
    /* pi = 4 * (a + b) ^ 2 / t */
    mpf_add(pi, a, b);
    mpf_mul(pi, pi, pi);
    mpf_mul_ui(t, t, 4);
    mpf_div(pi, pi, t);

    fprintf(stderr, "pi\n");

    mp_exp_t exponent;
    char *result = mpf_get_str(NULL, &exponent, 10, 0, pi);

    printf("%s\n", result + exponent);

    /* Clear all the mess :D */
    free(result);
    mpf_clear(pi);
    mpf_clear(a);
    mpf_clear(b);
    mpf_clear(t);
    mpf_clear(p);
    mpf_clear(a0);
    mpf_clear(b0);
    mpf_clear(t0);
    mpf_clear(p0);

    return 0;
}

